home *** CD-ROM | disk | FTP | other *** search
/ Tech Arsenal 1 / Tech Arsenal (Arsenal Computer).ISO / tek-19 / iritsm3s.zip / CBSP_AUX.C < prev    next >
C/C++ Source or Header  |  1991-06-02  |  10KB  |  276 lines
  1. /******************************************************************************
  2. * CBsp-Aux.c - Bspline curve auxilary routines.                      *
  3. *******************************************************************************
  4. * Written by Gershon Elber, Aug. 90.                          *
  5. ******************************************************************************/
  6.  
  7. #ifdef __MSDOS__
  8. #include <stdlib.h>
  9. #endif /* __MSDOS__ */
  10.  
  11. #include <ctype.h>
  12. #include <stdio.h>
  13. #include <string.h>
  14. #include "cagd_loc.h"
  15.  
  16. /******************************************************************************
  17. * Given a bspline curve - subdivide it into two at the given parametric value.*
  18. * Returns pointer to first curve in a list of two curves (subdivided ones).   *
  19. * The subdivision is achieved by inserting (order-1) knot at the given param. *
  20. * value t and spliting the control polygon and knot vector at that point.     *
  21. ******************************************************************************/
  22. CagdCrvStruct *BspCrvSubdivAtParam(CagdCrvStruct *Crv, CagdRType t)
  23. {
  24.     CagdBType
  25.     IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv);
  26.     int j,
  27.     k = Crv -> Order,
  28.     Len = Crv -> Length,
  29.         KVLen = k + Len,
  30.     Index1 = BspKnotLastIndexL(Crv -> KnotVector, KVLen, t),
  31.     Index2 = BspKnotFirstIndexG(Crv -> KnotVector, KVLen, t),
  32.     MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
  33.     CagdCrvStruct
  34.     *RefinedCrv = BspCrvKnotInsertNSame(Crv, t, k - 1),
  35.     *LCrv = BspCrvNew(Index1 + 1, k, Crv -> PType),
  36.     *RCrv = BspCrvNew(Len - Index2 + k, k, Crv -> PType);
  37.  
  38.     /* Copy Curve into LCrv. */
  39.     for (j = IsNotRational; j <= MaxCoord; j++)
  40.     GEN_COPY(LCrv -> Points[j],
  41.          RefinedCrv -> Points[j],
  42.                  sizeof(CagdRType) * (Index1 + 1));
  43.     GEN_COPY(LCrv -> KnotVector,
  44.          RefinedCrv -> KnotVector,
  45.          sizeof(CagdRType) * (Index1 + k));
  46.     /* Close the knot vector with multiplicity Order: */
  47.     LCrv -> KnotVector[Index1 + k] = LCrv -> KnotVector[Index1 + k - 1];
  48.  
  49.     /* Copy Curve into RCrv. */
  50.     for (j = IsNotRational; j <= MaxCoord; j++)
  51.     GEN_COPY(RCrv -> Points[j],
  52.          &RefinedCrv -> Points[j][Index1],
  53.          sizeof(CagdRType) * (Len - Index2 + k));
  54.     GEN_COPY(&RCrv -> KnotVector[1],
  55.          &RefinedCrv -> KnotVector[Index1 + 1],
  56.          sizeof(CagdRType) * (Len - Index2 + k + k - 1));
  57.     /* Make sure knot vector starts with multiplicity Order: */
  58.     RCrv -> KnotVector[0] = RCrv -> KnotVector[1];
  59.  
  60.     LCrv -> Pnext = RCrv;
  61.  
  62.     CagdCrvFree(RefinedCrv);
  63.  
  64.     return LCrv;
  65. }
  66.  
  67. /******************************************************************************
  68. * Return a new curve, identical to the original but with one degree higher    *
  69. ******************************************************************************/
  70. CagdCrvStruct *BspCrvDegreeRaise(CagdCrvStruct *Crv)
  71. {
  72.     CagdBType
  73.     IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv);
  74.     int i, i2, j, RaisedLen,
  75.     Order = Crv -> Order,
  76.     Length = Crv -> Length,
  77.     MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
  78.     CagdCrvStruct *RaisedCrv;
  79.  
  80.     if (Order > 2) {
  81.     FATAL_ERROR(CAGD_ERR_NOT_IMPLEMENTED);
  82.     return NULL;
  83.     }
  84.  
  85.     /* If curve is linear, degree raising means basically to increase the     */
  86.     /* knot multiplicity of each segment by one and add a middle point for    */
  87.     /* each such segment.                              */
  88.     RaisedLen = Length * 2 - 1;
  89.     RaisedCrv = BspCrvNew(RaisedLen, Order + 1, Crv -> PType);
  90.  
  91.     /* Update the control polygon; */
  92.     for (j = IsNotRational; j <= MaxCoord; j++)              /* First point. */
  93.     RaisedCrv -> Points[j][0] = Crv -> Points[j][0];
  94.  
  95.     for (i = 1, i2 = 1; i < Length; i++, i2 += 2)
  96.     for (j = IsNotRational; j <= MaxCoord; j++) {
  97.         RaisedCrv -> Points[j][i2] =
  98.         Crv -> Points[j][i-1] * 0.5 + Crv -> Points[j][i] * 0.5;
  99.         RaisedCrv -> Points[j][i2 + 1] = Crv -> Points[j][i];
  100.     }
  101.  
  102.     /* Update the knot vector. */
  103.     for (i = 0; i < 3; i++)
  104.     RaisedCrv -> KnotVector[i] = Crv -> KnotVector[0];
  105.  
  106.     for (i = 2, j = 3; i < Length; i++, j += 2)
  107.     RaisedCrv -> KnotVector[j] = RaisedCrv -> KnotVector[j + 1] = 
  108.         Crv -> KnotVector[i];
  109.  
  110.     for (i = j; i < j + 3; i++)
  111.     RaisedCrv -> KnotVector[i] = Crv -> KnotVector[Length];
  112.  
  113.     return RaisedCrv;
  114. }
  115.  
  116. /******************************************************************************
  117. * Return a normalized vector, equal to the tangent to Crv at parameter t.     *
  118. * Algorithm: insert (order - 1) knots and return control polygon tangent.     *
  119. ******************************************************************************/
  120. CagdVecStruct *BspCrvTangent(CagdCrvStruct *Crv, CagdRType t)
  121. {
  122.     static CagdVecStruct P2;
  123.     CagdVecStruct P1;
  124.     CagdRType TMin, TMax;
  125.     int k = Crv -> Order,
  126.     Len = Crv -> Length,
  127.     Index = BspKnotLastIndexL(Crv -> KnotVector, k + Len, t);
  128.     CagdPointType
  129.     PType = Crv -> PType;
  130.     CagdCrvStruct *RefinedCrv;
  131.  
  132.     if (!BspKnotParamInDomain(Crv -> KnotVector, Len, k, t))
  133.     FATAL_ERROR(CAGD_ERR_T_NOT_IN_CRV);
  134.  
  135.     BspCrvDomain(Crv, &TMin, &TMax);
  136.  
  137.     if (APX_EQ(t, TMin)) {
  138.     /* Use Crv starting tangent direction. */
  139.     CagdCoerceToE3(P1.Vec, Crv -> Points, 0, PType);
  140.     CagdCoerceToE3(P2.Vec, Crv -> Points, 1, PType);
  141.     }
  142.     else if (APX_EQ(t, TMax)) {
  143.     /* Use Crv ending tangent direction. */
  144.     CagdCoerceToE3(P1.Vec, Crv -> Points, Len - 2, PType);
  145.     CagdCoerceToE3(P2.Vec, Crv -> Points, Len - 1, PType);
  146.     }
  147.     else {
  148.     RefinedCrv = BspCrvKnotInsertNSame(Crv, t, k - 1);
  149.  
  150.     CagdCoerceToE3(P1.Vec, RefinedCrv -> Points, Index, PType);
  151.     CagdCoerceToE3(P2.Vec, RefinedCrv -> Points, Index + 1, PType);
  152.  
  153.     CagdCrvFree(RefinedCrv);
  154.     }
  155.  
  156.     CAGD_SUB_VECTOR(P2, P1);
  157.  
  158.     CAGD_NORMALIZE_VECTOR(P2);                /* Normalize the vector. */
  159.  
  160.     return &P2;
  161. }
  162.  
  163. /******************************************************************************
  164. * Return a normalized vector, equal to the binormal to Crv at parameter t.    *
  165. * Algorithm: insert (order - 1) knots and using 3 consecutive control points  *
  166. * at the refined location (p1, p2, p3), compute to binormal to be the cross   *
  167. * product of the two vectors (p1 - p2) and (p2 - p3).                  *
  168. ******************************************************************************/
  169. CagdVecStruct *BspCrvBiNormal(CagdCrvStruct *Crv, CagdRType t)
  170. {
  171.     static CagdVecStruct P3;
  172.     CagdVecStruct P1, P2;
  173.     CagdRType TMin, TMax;
  174.     int k = Crv -> Order,
  175.     Len = Crv -> Length,
  176.     Index = BspKnotLastIndexL(Crv -> KnotVector, k + Len, t);
  177.     CagdPointType
  178.     PType = Crv -> PType;
  179.     CagdCrvStruct *RefinedCrv;
  180.  
  181.     if (!BspKnotParamInDomain(Crv -> KnotVector, Len, k, t))
  182.     FATAL_ERROR(CAGD_ERR_T_NOT_IN_CRV);
  183.  
  184.     /* Can not compute for linear curves. */
  185.     if (k <= 2) return NULL;
  186.  
  187.     BspCrvDomain(Crv, &TMin, &TMax);
  188.  
  189.     if (APX_EQ(t, TMin)) {
  190.     /* Use Crv starting tangent direction. */
  191.     CagdCoerceToE3(P1.Vec, Crv -> Points, 0, PType);
  192.     CagdCoerceToE3(P2.Vec, Crv -> Points, 1, PType);
  193.     CagdCoerceToE3(P3.Vec, Crv -> Points, 2, PType);
  194.     }
  195.     else if (APX_EQ(t, TMax)) {
  196.     /* Use Crv ending tangent direction. */
  197.     CagdCoerceToE3(P1.Vec, Crv -> Points, Len - 3, PType);
  198.     CagdCoerceToE3(P2.Vec, Crv -> Points, Len - 2, PType);
  199.     CagdCoerceToE3(P3.Vec, Crv -> Points, Len - 1, PType);
  200.     }
  201.     else {
  202.     RefinedCrv = BspCrvKnotInsertNSame(Crv, t, k - 1);
  203.  
  204.     CagdCoerceToE3(P1.Vec, RefinedCrv -> Points, Index, PType);
  205.     CagdCoerceToE3(P2.Vec, RefinedCrv -> Points, Index + 1, PType);
  206.     CagdCoerceToE3(P3.Vec, RefinedCrv -> Points, Index + 2, PType);
  207.  
  208.     CagdCrvFree(RefinedCrv);
  209.     }
  210.  
  211.     CAGD_SUB_VECTOR(P1, P2);
  212.     CAGD_SUB_VECTOR(P2, P3);
  213.  
  214.     CROSS_PROD(P3.Vec, P1.Vec, P2.Vec);
  215.  
  216.     if ((t = CAGD_LEN_VECTOR(P3)) < EPSILON)
  217.     return NULL;
  218.     else
  219.     CAGD_DIV_VECTOR(P3, t);                /* Normalize the vector. */
  220.  
  221.     return &P3;
  222. }
  223.  
  224. /******************************************************************************
  225. * Return a normalized vector, equal to the normal to Crv at parameter t.      *
  226. * Algorithm: returns the cross product of the curve tangent and binormal.     *
  227. ******************************************************************************/
  228. CagdVecStruct *BspCrvNormal(CagdCrvStruct *Crv, CagdRType t)
  229. {
  230.     static CagdVecStruct N, *T, *B;
  231.  
  232.     T = BspCrvTangent(Crv, t);
  233.     B = BspCrvBiNormal(Crv, t);
  234.  
  235.     if (T == NULL || B == NULL) return NULL;
  236.  
  237.     CROSS_PROD(N.Vec, T -> Vec, B -> Vec);
  238.  
  239.     CAGD_NORMALIZE_VECTOR(N);                /* Normalize the vector. */
  240.  
  241.     return &N;
  242. }
  243.  
  244. /******************************************************************************
  245. * Return a new curve, equal to the derived curve. If the original curve is    *
  246. * rational, NULL is return (curve must be non-rational for this one).          *
  247. * Let old control polygon be P(i), i = 0 to k-1, and Q(i) be new one then:    *
  248. * Q(i) = (k - 1) * P(i+1) - P(i), i = 0 to k-2.                      *
  249. ******************************************************************************/
  250. CagdCrvStruct *BspCrvDerive(CagdCrvStruct *Crv)
  251. {
  252.     CagdBType
  253.     IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv);
  254.     int i, j,
  255.     k = Crv -> Order,
  256.     Len = Crv -> Length,
  257.     MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
  258.     CagdCrvStruct *DerivedCrv;
  259.  
  260.     if (!IsNotRational || k < 3)
  261.     FATAL_ERROR(CAGD_ERR_RAT_LIN_NO_SUPPORT);
  262.  
  263.     DerivedCrv = BspCrvNew(Len - 1, k - 1, Crv -> PType);
  264.  
  265.     for (i = 0; i < Len - 1; i++)
  266.     for (j = IsNotRational; j <= MaxCoord; j++)
  267.         DerivedCrv -> Points[j][i] =
  268.         (k - 1) * (Crv -> Points[j][i+1] - Crv -> Points[j][i]);
  269.  
  270.     GEN_COPY(DerivedCrv -> KnotVector,
  271.              &Crv -> KnotVector[1],
  272.              sizeof(CagdRType) * (k + Len - 2));
  273.  
  274.     return DerivedCrv;
  275. }
  276.